This script shows results of our analysis of NDWI trends
(insignificant trends) and congruency with ERA5 climate trend
patterns.
### ------------------------------------------------------------------------------------------------------------ ###
# STEP 1: function, libraries and plot preparations
### ------------------------------------------------------------------------------------------------------------ ###
#define function (cite: CRAN 2014: USGS Sen's Slope: https://rdrr.io/github/USGS-R/smwrStats/man/senSlope.htm)
slope_fun=function(r){
senSlope <- function(formula, data, subset, na.action, intercept='Ac',
CI=.95) {
# Coding History:
# 2000Oct27 JRSlack Original coding as kensen
# 2011May02 DLLorenz Conversion to R--as sen slope only
# 2011Oct25 DLLorenz Update for package
# 2013Apr30 DLLorenz Bug fixes
# 2014Dec29 DLLorenz Conversion to roxygen header
##
## Define the variance of Kendall's S function, required for conf. int.
vark <- function(y, x) {
ties.y <- rle(sort(y))$lengths
ties.x <- rle(sort(x))$lengths
n <- length(y)
t1 <- n * (n - 1) * (2 * n + 5)
ty2 <- sum(ties.y * (ties.y - 1) * (2 * ties.y + 5))
tx2 <- sum(ties.x * (ties.x - 1) * (2 * ties.x + 5))
v <- (t1 - ty2 -tx2)/18
v1 <- sum(ties.y * (ties.y - 1)) * sum(ties.x * (ties.x - 1)) /
(2 * n * (n - 1))
v2 <- sum(ties.y * (ties.y - 1) * (ties.y - 2)) *
sum(ties.x * (ties.x - 1) * (ties.x - 2)) /
(9 * n * (n - 1) * (n - 2))
v <- v + v1 + v2
return (v)
}
## Process formula as with regular linear regression!
call <- match.call()
m <- match.call(expand.dots = FALSE)
m$intercept <- m$CI <- NULL
m$drop.unused.levels <- TRUE
m[[1]] <- as.name("model.frame")
m <- eval(m, sys.parent())
if(ncol(m) != 2L)
stop("Must specify a single response and a single explanatory variable in the formula")
## Extract missing value info
na.action <- attr(m, "na.action")
y <- m[, 1L]
yname <- names(m)[1L]
x <- m[,2]
xname <- names(m)[2L]
n <- length(y)
if (n < 3L)
stop("model data should effectively be longer than 2.")
## Calculate the statistics
## A fast, efficient way to compute all possible slopes:
slopes <- unlist(lapply(seq(along = y), function(i, y, t)
((y[i] - y[1L:i]) / (t[i] - t[1L:i])), y, x))
slopes <- sort(slopes) # removes missings (due to ties in x)
sen.slope <- median(slopes)
## Compute the variance of S, accounting only for ties in x
varS <- vark(seq(along=y), x) # Forces no ties in y
Calpha <- -qnorm((1-CI)/2) * sqrt(varS)
M1 <- as.integer((length(slopes) - Calpha)/2)
M2 <- as.integer((length(slopes) + Calpha)/2) + 1
sen.CI <- slopes[c(M1, M2)]
names(sen.CI) <- paste(c("Lower", "Upper"),
substring(sprintf("%.2f", CI), 2), sep="") # drop 0
## Median of the data values.
median.y <- median(y)
## Median of the time values.
median.x <- median(x)
## A line representing the trend of the data then is given by
##
## y(t) = sen.slope*(t-med.time)+med.data
##
## This line has the slope of the trend and passes through
## the point (t=med.time, y=med.data)
## Compute the coefficients for the line: intercept and slope
coef <- c(sen.slope*(-median.x)+median.y, sen.slope)
fits <- coef[1L] + coef[2L]*x
resids <- y - fits
if(intercept == 'A1m') {
coef[1L] <- coef[1L] + median(resids)
resids <- y - coef[1L] - coef[2L]*x
}
## Return the statistics.
retval <- list(call=call, coefficients=coef, slope.CI=sen.CI, residuals=resids,
fitted.values=fits, na.action=na.action, x=x, y=y, var.names=c(yname, xname),
model=m)
oldClass(retval) <- "senSlope"
return(retval)
}
ts_cell=as.data.frame(r)
df_ts_cell=data.frame(cbind('year'=c(seq(1982, 2022)),ts_cell))
#print(df_ts_cell)
if(length(which(is.na(df_ts_cell[,2])))<40){
Sslope_vsm_USGS=senSlope(df_ts_cell[,2] ~ df_ts_cell[,1], data=df_ts_cell, na.action = na.exclude)
return(Sslope_vsm_USGS$coefficients[2])
}else{
return(NA)
}
}
#load libraries
library(terra)
library(sf)
library(tidyr)
library(dplyr)
library(stringr)
library(ggplot2)
library(ggspatial)
library(tidyterra)
library(ggpubr)
library(gridExtra)
library(plotly)
library(rnaturalearth)
library(rnaturalearthdata)
#theme for map plots
myTheme <- theme(legend.text = element_text(size = 10),
legend.title = element_text(size = 10),
legend.key.size = unit(1, 'cm'),
legend.position = 'bottom',
panel.grid.major = element_blank(),
panel.background = element_rect(fill = "#cccccc"))
pal <- colorRampPalette(c('goldenrod3','goldenrod3','goldenrod3',"#f7c973","#f7c973", "#f7bf59", '#ecebeb', 'lightblue', "#209fb6", "#209fb6","#209fb6",'darkblue','darkblue','darkblue'))
pal_t= colorRampPalette(c('blue', 'darkblue','white','#ffd000','#ff5100'))
##create base map with countries of europe - spatial polygons
#basemap
world_map <- rnaturalearth::ne_countries(scale = 50, returnclass = 'sf')
european_union <- c("Austria","Albania","Andorra","Bosnia and Herz.", "Belgium","Bulgaria","Croatia","Cyprus",
"Czech Rep.","Denmark","Estonia","Finland","France",
"Germany","Greece","Hungary","Ireland","Italy","Kosovo","Latvia",
"Lithuania","Luxembourg","Malta","Montenegro", "Netherlands","Poland",
"Portugal","Romania","Slovakia","Slovenia","Spain", "Serbia",
"Sweden","United Kingdom","Norway", "Iceland", "Switzerland", "Czechia", "Macedonia")
european_union_map <-
world_map %>%
filter(name %in% european_union)
bbox_europe <- st_bbox(c(xmin = -30, ymin = 30, xmax = 35, ymax = 70), crs = st_crs(4326))
#european_union_map_proj=st_transform(european_union_map, crs=st_crs(bbox_europe))
european_union_map_cropped <- st_crop(european_union_map, bbox_europe)
df <-
tibble(country = european_union,
some_value = runif(length(european_union)))
map <-
european_union_map_cropped %>%
left_join(df, by = c("name" = "country"))
STEP 2: Calculate Sen’s Slope for ERA5 climate data
#load climate time series ERA5 (temperature & precipitation in one data set)
VSM_rast=rast('/home/laurag/Arbeit/wwu/R/workspaces/satellite_ts_R/random_sample_analysis/copernicus/PuT/data.nc')
#check which layers represent temperature (here: 1:492) and precipitation (here: 493:984)
#time(VSM_rast)
#names(VSM_rast)
#temperature
era5_var_temp=VSM_rast[[1:492]]
#precipitation
era5_var_precip=VSM_rast[[493:984]]
#generate variable for monthly filtering
patterns=c(sapply(str_pad(seq(12), 2, pad = "0"), function(s){paste0('-',s,'-01')}))
#select Jun-Aug
patterns=c(patterns[6],patterns[7],patterns[8])
### ------------------------------------------------------------------------------------------------------------ ###
#calculate and plot temperature trends
### ------------------------------------------------------------------------------------------------------------ ###
#app() trend function, use multiple cores
era5_trend_tif_list_temp=lapply(patterns, function(p){
era5_mon_temp=era5_var_temp[[grep(pattern=p,time(era5_var_temp))]]
test_app_out=terra::app(era5_mon_temp,slope_fun, cores =4)
#write climate trend raster
#terra::writeRaster(test_app_out, paste0('/home/laurag/Arbeit/wwu/R/workspaces/satellite_ts_R/random_sample_analysis/copernicus/PuT/era5_trend_tp_m_',substr(p, start=2, stop=3),'.tif'))
return(test_app_out)
})
season_r_temp=rast(era5_trend_tif_list_temp)
e=ext(-25, 34, 46, 72)
season_r_temp_cr=crop(season_r_temp, e)
season_r_temp_mean=app(season_r_temp_cr, mean)
names(season_r_temp_mean)=c('Summer (JJA mean)')
ggplot() +
geom_spatraster(data=season_r_temp_mean) +
geom_sf(data = map, color = "grey30", fill = NA, size=0.5) +
facet_wrap(~lyr, ncol=2) +
coord_sf(xlim = c(-25, 30), ylim = c(47, 70))+
myTheme+ theme(panel.grid.major = element_blank(), panel.background = element_rect(fill = "#aca4a4", colour = 'grey'))+
scale_fill_gradientn(colours = pal_t(7), limits= c(-0.22,0.22), breaks =c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15), na.value = "transparent") +
labs(fill = "Temperature Change [°C]/yr")+
ggtitle('Temperature Trends')
Summer (JJA mean)
20
°
W
10
°
W
0
°
10
°
E
20
°
E
30
°
E
50
°
N
55
°
N
60
°
N
65
°
N
70
°
N
-0.15
-0.10
-0.05
0.00
0.05
0.10
0.15
Temperature Change [°C]/yr
Temperature Trends
### ------------------------------------------------------------------------------------------------------------ ###
#calculate and plot precipitation trends
### ------------------------------------------------------------------------------------------------------------ ###
era5_trend_tif_list_precip=lapply(patterns, function(p){
era5_mon_precip=era5_var_precip[[grep(pattern=p,time(era5_var_precip))]]
test_app_out=terra::app(era5_mon_precip,slope_fun, cores =4)
#write climate trend raster
#terra::writeRaster(test_app_out, paste0('/home/laurag/Arbeit/wwu/R/workspaces/satellite_ts_R/random_sample_analysis/copernicus/PuT/era5_trend_tp_m_',substr(p, start=2, stop=3),'.tif'))
return(test_app_out)
})
season_r_precip=rast(era5_trend_tif_list_precip)
e=ext(-25, 34, 46, 72)
season_r_precip_cr=crop(season_r_precip, e)
season_r_precip_mean=app(season_r_precip_cr, mean)
names(season_r_precip_mean)=c('Summer (JJA mean)')
# convert unit from m to mm
season_r_precip_mean_mm=season_r_precip_mean*1000
ggplot() +
#coord_sf(xlim = c(-30, 33), ylim = c(37, 73))+
geom_spatraster(data=season_r_precip_mean_mm) +
geom_sf(data = map, color = "grey30", fill = NA, size=0.5) +
facet_wrap(~lyr, ncol=2) +
coord_sf(xlim = c(-25, 30), ylim = c(47, 70))+
myTheme+ theme(panel.grid.major = element_blank(), panel.background = element_rect(fill = "#aca4a4", colour = 'grey'))+
scale_fill_gradient2(low='goldenrod', mid= 'white', high='darkblue',midpoint= 0, limits= c(-0.06,0.06), breaks =c(-0.05, -0.025, 0, 0.025, 0.05), na.value = "transparent") +
labs(fill = "Tot. Precipitation Change [mm]/yr")+
ggtitle('Precipitation Trends')
Summer (JJA mean)
20
°
W
10
°
W
0
°
10
°
E
20
°
E
30
°
E
50
°
N
55
°
N
60
°
N
65
°
N
70
°
N
-0.050
-0.025
0.000
0.025
0.050
Tot. Precipitation Change [mm]/yr
Precipitation Trends
#ggsave(file = '/home/laurag/Arbeit/wwu/R/workspaces/satellite_ts_R/random_sample_analysis/copernicus/PuT/EU_tp_m_ERA5_SWmean.png', device='png',
# width = 14, # The width of the plot in inches
# height = 10)
STEP 3: Prepare NDWI trend data set for comparision with climate
trends
(also includeded in script ‘mtEUpeat_significantTrends_p2’)
### ------------------------------------------------------------------------------------------------------------ ###
#load data and filter significant trends
### ------------------------------------------------------------------------------------------------------------ ###
#define input path
path_trend_points='/home/laurag/Arbeit/wwu/R/workspaces/satellite_ts_R/random_sample_analysis/output/'
#load trend points
trendstats_filtered_outliers=readRDS(paste0(path_trend_points, 'trendpoints_fwo_filtered_list.rds'))
#filter NDWI trend dataset on count of timestamps per time-series > 10
#and p-value <0.05 (95% confidence level) for May-October
slope_filtered_sig=lapply(seq(5,10), function(i){
m=month.name[i]
points_filtered= trendstats_filtered_outliers[[i]] %>% filter(.[[1]] > 10) %>% filter(.[[3]]<= 0.05)
trendstats_sig=as.data.frame(cbind(c(rep(m, nrow(points_filtered))),points_filtered))
colnames(trendstats_sig)=c('mon', 'count', 'tau', 'p_val', 'Sint', 'Sslope', 'geom')
trendstats_sig_sf=st_as_sf(trendstats_sig, geom=trendstats_sig$geom)
return(trendstats_sig_sf)
})
#combine in one simple feature collection (long format)
slope_filtered_sig_bind=do.call(rbind,slope_filtered_sig)
slope_filtered_sig_bind$mon=factor(slope_filtered_sig_bind$mon, levels=month.name[5:10])
### ------------------------------------------------------------------------------------------------------------ ###
##calculate mean trend of summer months mean
### ------------------------------------------------------------------------------------------------------------ ###
#extract common coordinates (Jun-Aug) for points with significant trends
equal_points67=slope_filtered_sig[[2]][which(st_geometry(slope_filtered_sig[[2]]) %in% st_geometry(slope_filtered_sig[[3]])),]
equal_points672=equal_points67[which(st_geometry(equal_points67) %in% st_geometry(slope_filtered_sig[[4]])),]
Sslope_jun_common=slope_filtered_sig[[2]][which(st_geometry(slope_filtered_sig[[2]]) %in% st_geometry(equal_points672)),]
Sslope_jul_common=slope_filtered_sig[[3]][which(st_geometry(slope_filtered_sig[[3]]) %in% st_geometry(equal_points672)),]
Sslope_aug_common=slope_filtered_sig[[4]][which(st_geometry(slope_filtered_sig[[4]]) %in% st_geometry(equal_points672)),]
summer_months_bind=cbind(Sslope_jul_common,Sslope_jun_common,Sslope_aug_common)
summer_months_bind_df=as.data.frame(cbind(as.data.frame(summer_months_bind[,6])[,1], as.data.frame(summer_months_bind[,12])[,1], as.data.frame(summer_months_bind[,18])))
colnames(summer_months_bind_df)=c('Jun','Jul','Aug','geometry')
summer_months_bind_df[,1]=as.numeric(summer_months_bind_df[,1])
summer_months_bind_df[,2]=as.numeric(summer_months_bind_df[,2])
summer_months_bind_df[,3]=as.numeric(summer_months_bind_df[,3])
#calculate mean
summer_months_mean = lapply(seq(1,nrow(summer_months_bind_df)), function(r){
mean_ft=mean(c(summer_months_bind_df$Jul[r],summer_months_bind_df$Jun[r],summer_months_bind_df$Aug[r]), na.rm=T)
return(mean_ft)
})
summer_months_mean_bind=as.data.frame(do.call(rbind,summer_months_mean))
summer_months_mean_cbind=cbind(summer_months_mean_bind, summer_months_bind_df[,4])
summer_months_mean_cbind_df=as.data.frame(summer_months_mean_cbind)
summer_months_mean_cbind_sf=st_as_sf(summer_months_mean_cbind_df, geom=summer_months_mean_cbind_df$geometry)
summer_months_mean_cbind_sf=summer_months_mean_cbind_sf[,-2]
#plot equal points summer months
summer_months_mean_cbind_sf$V1=as.numeric(summer_months_mean_cbind_sf$V1)
ggplot() +
geom_sf(data=map, color='darkgrey',fill = 'white', size=0.05) +
geom_sf(data =summer_months_mean_cbind_sf, aes(colour=summer_months_mean_cbind_sf$V1), size = .1) +
#facet_wrap(~Month)+
coord_sf(xlim = c(-30, 33), ylim = c(47, 71)) +
scale_colour_stepsn(colours = pal(7), limits= c(-0.006,0.006), breaks =c(-0.006,-0.003,-0.001,0,0.001,0.003,0.006))+
labs(colour = "SensSlope") +
myTheme+
ggtitle('Summer (JJA) mean trend')
50
°
N
55
°
N
60
°
N
65
°
N
70
°
N
30
°
W
20
°
W
10
°
W
0
°
10
°
E
20
°
E
30
°
E
-0.006
-0.003
-0.001
0.000
0.001
0.003
0.006
SensSlope
Summer (JJA) mean trend
STEP 4: Congruence Analysis of NDWI trends and climatic trends
###------------------------------------------------------------------------------------------------------###
### Test for concruence of NDWI and precipitation trends
###------------------------------------------------------------------------------------------------------###
#summary(summer_months_mean_cbind_sf)
#extract clim trend at point location
trendpoints_pixel_extracted_precip=terra::extract(season_r_precip_mean_mm, summer_months_mean_cbind_sf, bind = TRUE)
trendpoints_pixel_extracted_precip_sf=st_as_sf(trendpoints_pixel_extracted_precip)
#p=5
match_var_all_precip=lapply(seq(nrow(trendpoints_pixel_extracted_precip)), function(p){
if(isTRUE(as.numeric(trendpoints_pixel_extracted_precip_sf[p,1])[1] > 0) & isTRUE(as.numeric(trendpoints_pixel_extracted_precip_sf[p,2])[1]>0)){
match_var = 'pos'
}else if(isTRUE(as.numeric(trendpoints_pixel_extracted_precip_sf[p,1])[1]<0) & isTRUE(as.numeric(trendpoints_pixel_extracted_precip_sf[p,2])[1]<0)){
match_var = 'neg'
}else{
match_var = 'none'
}
return(match_var)
})
match_var_all_precip_bind=as.data.frame(do.call(rbind,match_var_all_precip))
trendpoints_pixel_extracted_precip_matched=cbind(trendpoints_pixel_extracted_precip_sf,match_var_all_precip_bind)
colnames(trendpoints_pixel_extracted_precip_matched)=c('Sslope_NDWI','Sslope_precip', 'matched', 'geometry')
trendpoints_pixel_extracted_matched_p=trendpoints_pixel_extracted_precip_matched
#st_write(trendpoints_pixel_extracted_precip_matched, '/home/laurag/Arbeit/wwu/R/workspaces/satellite_ts_R/random_sample_analysis/copernicus/SslopePrecip_SslopeNDWI_match.gpkg')
trendpoints_pixel_extracted_matched_na_rem_p = as.data.frame(trendpoints_pixel_extracted_matched_p)[complete.cases(as.data.frame(trendpoints_pixel_extracted_matched_p)[,1:3]),]
colnames(trendpoints_pixel_extracted_matched_na_rem_p)=c('Sslope_NDWI','Sslope_Prec', 'matched_Prec', 'geometry')
###------------------------------------------------------------------------------------------------------###
### Test for concruence of NDWI and temperature trends
###------------------------------------------------------------------------------------------------------###
trendpoints_pixel_extracted_temp=terra::extract(season_r_temp_mean, summer_months_mean_cbind_sf, bind = TRUE)
trendpoints_pixel_extracted_temp_sf=st_as_sf(trendpoints_pixel_extracted_temp)
match_var_all_temp=lapply(seq(nrow(trendpoints_pixel_extracted_temp_sf)), function(p){
if(isTRUE(as.numeric(trendpoints_pixel_extracted_temp_sf[p,1])[1] > 0) & isTRUE(as.numeric(trendpoints_pixel_extracted_temp_sf[p,2])[1]<0.05)){
match_var = 'pos'
}else if(isTRUE(as.numeric(trendpoints_pixel_extracted_temp_sf[p,1])[1]<0) & isTRUE(as.numeric(trendpoints_pixel_extracted_temp_sf[p,2])[1]>0.05)){
match_var = 'neg'
}else{
match_var = 'none'
}
return(match_var)
})
match_var_all_temp_bind=as.data.frame(do.call(rbind,match_var_all_temp))
trendpoints_pixel_extracted_temp_matched=cbind(trendpoints_pixel_extracted_temp_sf,match_var_all_temp_bind)
colnames(trendpoints_pixel_extracted_temp_matched)=c('Sslope_NDWI','Sslope_temp', 'matched', 'geometry')
trendpoints_pixel_extracted_matched_t=trendpoints_pixel_extracted_temp_matched
trendpoints_pixel_extracted_matched_na_rem_t = as.data.frame(trendpoints_pixel_extracted_matched_t)[complete.cases(as.data.frame(trendpoints_pixel_extracted_matched_t)[,1:3]),]
colnames(trendpoints_pixel_extracted_matched_na_rem_t)=c('Sslope_NDWI','Sslope_Temp', 'matched_Temp', 'geometry')
###---------------------------------------------------------------------------------------------------------###
### Combine results of congruence analysis between NDWI trends, temperature and precipitation to one data frame
###---------------------------------------------------------------------------------------------------------###
tot_lm_ndwiclim=cbind(trendpoints_pixel_extracted_matched_na_rem_t,trendpoints_pixel_extracted_matched_na_rem_p[,2:3])
colnames(tot_lm_ndwiclim)## [1] "Sslope_NDWI" "Sslope_Temp" "matched_Temp" "geometry" "Sslope_Prec"
## [6] "matched_Prec"#colnames(tot_lm_ndwiclim_res)=c('Sslope_NDWI','Sslope_Temp', 'matched_Temp', 'resid_temp', 'Sslope_Prec', 'matched_Prec', 'resid_Prec', 'geometry', 'geom')
tot_lm_ndwiclim_sf=st_as_sf(tot_lm_ndwiclim, geom = tot_lm_ndwiclim$geometry)
#preapare spatial feature data set (congruence analysis) for plots
pos_neg_congr=lapply(seq(nrow(tot_lm_ndwiclim_sf)), function(r){
if(as.numeric(tot_lm_ndwiclim_sf$Sslope_NDWI[r])[1]>0){
congr_prec=paste0(tot_lm_ndwiclim_sf$matched_Prec[r],'_pos')
congr_temp=paste0(tot_lm_ndwiclim_sf$matched_Temp[r],'_pos')
return(c(congr_prec,congr_temp))
}else if(as.numeric(tot_lm_ndwiclim_sf$Sslope_NDWI[r])[1]<0){
congr_prec=paste0(tot_lm_ndwiclim_sf$matched_Prec[r],'_neg')
congr_temp=paste0(tot_lm_ndwiclim_sf$matched_Temp[r],'_neg')
return(c(congr_prec,congr_temp))
}else{
return(c(NA,NA))
}
})
congr_vec_bind=do.call(rbind,pos_neg_congr)
tot_lm_ndwiclim_pn=cbind(tot_lm_ndwiclim_sf,congr_vec_bind)
tot_lm_ndwiclim_pn=tot_lm_ndwiclim_pn[,-9]
colnames(tot_lm_ndwiclim_pn)=c(colnames(tot_lm_ndwiclim_pn)[1:5], 'congr_p', 'congr_t', 'geometry')
#colour experiments
#pal_error=colorRampPalette(rev(c('#ff3000','#fe7f00','#fe7f00','#fe7f00','#fe7f00', "#ffcb64", '#f8e1b3','#f8e1b3','#f8e1b3','#ffcb64', '#fe7f00', '#fe7f00','#fe7f00','#fe7f00','#ff3000')))
#pal_error=scale_colour_manual(values=c("goldenrod3", "#752404","#209fb6", "#3e1558"))
#seperate positive and negative NDWI trends
pos_filtered=tot_lm_ndwiclim_pn %>% filter(Sslope_NDWI>0)
neg_filtered=tot_lm_ndwiclim_pn %>% filter(Sslope_NDWI<0)
###---------------------------------------------------------------------------------------------------------###
### Generate Congruency Plots
###---------------------------------------------------------------------------------------------------------###
## plot temperature
ggplot() + myTheme +
geom_sf(data=map, color='black' ,fill = 'white', lwd=0.5) +
#geom_sf(data =plot_sf, aes(colour=plot_sf[[1]]), size =3) + #c(data.frame(slope_filtered[[1]][,1])[,1]))
geom_sf(data=pos_filtered, aes(colour=pos_filtered$congr_t), size =0.8) +
#geom_sf(data =plot_sf, fill = 'white', col='black') +
coord_sf(xlim = c(-22, 29), ylim = c(42, 71))+
#scale_colour_stepsn(colours = pal_error(7), limits= c(-0.02,0.02), breaks =c(-0.02,-0.01,-0.001,0,0.001,0.01,0.02)) +
scale_colour_manual(values=c("hotpink4","#209fb6"))+
labs(colour = "Congruence")+
ggtitle('Temperature')
45
°
N
50
°
N
55
°
N
60
°
N
65
°
N
70
°
N
20
°
W
10
°
W
0
°
10
°
E
20
°
E
30
°
E
Congruence
none_pos
pos_pos
Temperature
ggplot() + myTheme +
geom_sf(data=map, color='black' ,fill = 'white', lwd=0.5) +
#geom_sf(data =plot_sf, aes(colour=plot_sf[[1]]), size =3) + #c(data.frame(slope_filtered[[1]][,1])[,1]))
geom_sf(data=neg_filtered, aes(colour=neg_filtered$congr_t), size =0.8) +
#geom_sf(data =plot_sf, fill = 'white', col='black') +
coord_sf(xlim = c(-22, 29), ylim = c(42, 71))+
#scale_colour_stepsn(colours = pal_error(7), limits= c(-0.02,0.02), breaks =c(-0.02,-0.01,-0.001,0,0.001,0.01,0.02)) +
scale_colour_manual(values=c("goldenrod3","hotpink4"))+
labs(colour = "Congruence")+
ggtitle('Temperature')
45
°
N
50
°
N
55
°
N
60
°
N
65
°
N
70
°
N
20
°
W
10
°
W
0
°
10
°
E
20
°
E
30
°
E
Congruence
neg_neg
none_neg
Temperature
## plot precipitation
ggplot() + myTheme +
geom_sf(data=map, color='black' ,fill = 'white', lwd=0.5) +
#geom_sf(data =plot_sf, aes(colour=plot_sf[[1]]), size =3) + #c(data.frame(slope_filtered[[1]][,1])[,1]))
geom_sf(data=pos_filtered, aes(colour=pos_filtered$congr_p), size =0.8) +
#geom_sf(data =plot_sf, fill = 'white', col='black') +
coord_sf(xlim = c(-22, 29), ylim = c(42, 71))+
#scale_colour_stepsn(colours = pal_error(7), limits= c(-0.02,0.02), breaks =c(-0.02,-0.01,-0.001,0,0.001,0.01,0.02)) +
scale_colour_manual(values=c("hotpink4","#209fb6"))+
labs(colour = "Congruence")+
ggtitle('Precipitation')
45
°
N
50
°
N
55
°
N
60
°
N
65
°
N
70
°
N
20
°
W
10
°
W
0
°
10
°
E
20
°
E
30
°
E
Congruence
none_pos
pos_pos
Precipitation
ggplot() + myTheme +
geom_sf(data=map, color='black' ,fill = 'white', lwd=0.5) +
#geom_sf(data =plot_sf, aes(colour=plot_sf[[1]]), size =3) + #c(data.frame(slope_filtered[[1]][,1])[,1]))
geom_sf(data=neg_filtered, aes(colour=neg_filtered$congr_p), size =0.8) +
#geom_sf(data =plot_sf, fill = 'white', col='black') +
coord_sf(xlim = c(-22, 29), ylim = c(42, 71))+
#scale_colour_stepsn(colours = pal_error(7), limits= c(-0.02,0.02), breaks =c(-0.02,-0.01,-0.001,0,0.001,0.01,0.02)) +
scale_colour_manual(values=c("goldenrod3","hotpink4"))+
labs(colour = "Congruence")+
ggtitle('Precipitation')
45
°
N
50
°
N
55
°
N
60
°
N
65
°
N
70
°
N
20
°
W
10
°
W
0
°
10
°
E
20
°
E
30
°
E
Congruence
neg_neg
none_neg
Precipitation
#ggsave('/home/laurag/Arbeit/wwu/R/workspaces/satellite_ts_R/random_sample_analysis/output/new_mean_summer_NDWItemp_pos_congr.svg',
#width = 10, height = 10)